Planetesimal disk evolution driven by planetesimal-planetesimal gravitational 

scattering. 

R. R. Rafikov 

Princeton University Observatory, Princeton, NJ 08544 
rrr@astro . princeton . edu 

ABSTRACT 

We investigate the process of an inhomogeneous planetesimal disk evolution caused 
by the planetesimal-planetesimal gravitational scattering. We develop a rather general 
approach based on the kinetic theory which self-consistently describes the evolution in 
time and space of both the disk surface density and its kinematic properties — dis- 
persions of eccentricity and inclination. The gravitational scattering of planetesimals 
is assumed to be in the dispersion-dominated regime which considerably simplifies an- 
alytical treatment. The resultant equations are of advection-diffusion type. Distance 
dependent scattering coefficients entering these equations are calculated analytically 
under the assumption of two-body scattering in the leading order in Coulomb loga- 
rithm. They are essentially nonlocal in nature. Our approach allows one to explore 
the dynamics of nonuniform planetesimal disks with arbitrary mass and random veloc- 
ity distributions. It can also naturally include other physical mechanisms which are 
important for the evolution of such disks — gas drag, migration, and so on. 

Subject headings: planets and satellites: general — solar system: formation — (stars:) 
planetary systems 



1. Introduction. 

The discovery of extrasolar giant planets around nearby stars (Mayor & Queloz 1995; Marcy 
et al. 2000; Vogt et al. 2000; Butler et al. 2001) has been one of the most exciting astrophysical 
findings of the last decade. It has revived interest in planetary sciences and stimulated many 
theoretical studies. One of the areas which has received a lot of attention recently is the formation 
of terrestrial planets — planets like Earth or Venus. 

There are many reasons for this interest. First of all, our own Solar System hosts several 
terrestrial planets and we must understand their formation mechanisms if we want to know the 
history of our planetary system. Second, despite the huge differences in their physical properties, 
giant planets are probably linked to the terrestrial ones through their formation mechanism, since 
it is widely believed that giant planets form as a result of gas accretion on a preexisting rocky core 



- 2 - 



which essentially was a massive Earth-type planet (Mizuno 1980; Stevenson 1982; Bodenheimer & 
Pollack 1986; see however Boss et al. 2002). Thus terrestrial planets are probably an important 
component in the evolutionary history of the giant planets. Third, a lot of effort is currently being 
directed toward designing and building space missions with the ultimate goal of detecting Earth- 
type planets around other stars. Theoretical understanding of how terrestrial planets had come 
into being will help us plan these missions most effectively. Finally, observations of IR emission 
from the dust and debris disks around nearby stars (Heinrichsen et al. 1999; Schneider 2001) must 
be interpreted in the context of the formation and evolution of such disks. They are very likely 
to be the outcome of the same collisional fragmentation and accumulation of massive rocky or icy 
bodies that form terrestrial planets. 

It is widely believed that the formation of Earth-type planets proceeded via agglomeration of 
large number of planetesimals — asteroid or comet-like rocky or icy bodies. The theory of this 
process was pioneered by Safronov (1972) who was the first to point out the importance of the evo- 
lution of the dynamical properties of the planetesimal disk for the evolution of its mass distribution. 
The discovery of a rapid "runaway" mode of planetary accretion by Wetherill & Stewart (1989) has 
made this issue even more important since this mechanism relies on the redistribution of the energy 
of planetesimal epicyclic motion between populations with different masses (Wetherill & Stewart 
1993; Kenyon & Luu 1998; Inaba et al. 2001). Significant progress has been made recently in 
understanding the dynamical evolution of homogeneous planetesimal disks — disks where gradients 
of surface density or dynamical properties (such as the dispersions of eccentricity and inclination of 
various planetesimal populations) are absent (Ida 1990; Stewart &; Ida 2000, hereafter SI; Ohtsuki 
et al. 2002). This assumption should be very good during the initial stages of planetesimal growth 
when there are no massive bodies in the disk. However as coagulation proceeds and planetary em- 
bryos — precursors of terrestrial planets — emerge this assumption runs into problems. It was first 
demonstrated by Ida & Makino (1993) using N-body simulations and then confirmed using orbit 
integrations (Tanaka & Ida 1997, 1999) that under a variety of conditions massive protoplanetary 
embryos would tend to repel planetesimal orbits, clearing out an annular gap in the disk around 
them. This process introduces a new spatial dimension into the problem and makes it much harder 
to treat. 

N-body simulations are not good tools to study the details of this process. The primary 
reason is that they become too time-consuming when one needs to follow the spatial and dynamical 
properties of a many-body system for many orbital periods. Moreover, the number of planetesimals 
which they can handle is not very large (less than 10 4 ) which precludes the consideration of realistic 
planetesimal disks containing huge number of bodies with masses spanning an enormous range - 
from one meter rocks to 100 km planetesimals. Finally, one would need to perform a large number 
of such simulations to explore the whole space of physical parameters relevant for protoplanetary 
disk evolution. Similar problems, although less severe, plague the performance of methods based 
of the direct integration of planetesimal orbits using some simplifying assumptions (Tanaka & Ida 
1997, 1999). 
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The large number of bodies under consideration poses no problem for the methods of kinetic 
theory. In this approach the disk is split into a set of planetesimal populations and each of them is 
characterized by three functions of position: the surface density and dispersions of eccentricity and 
inclination. Application of this method for studying planetesimal coagulation and the evolution of 
dynamical properties in homogeneous disks has proven to be very useful (Wetherill & Stewart 1989; 
Kenyon & Luu 1998). It can also be easily extended to the case of inhomogeneous disks by allowing 
the planetesimal properties to vary within the disk. This approach allows one to study many 
physical mechanisms important for the coagulation process — gravitational interactions between 
planetesimals increasing their random motion, gas drag and inelastic collisions damping them, 
migration, fragmentation, etc. The use of this statistical approach for exploring inhomogeneous 
disks was first undertaken by Petit & Henon (1987a, 1987b, 1988) in their studies of planetary 
rings. Physical conditions in planetary rings (high optical depth and frequent inelastic collisions) 
are very different from those which are thought to exist in protoplanetary disks. Thus we cannot 
directly apply the methods of Petit & Henon to study planetesimal disks, but the spirit in which 
their investigation was carried out can be largely preserved. 

In this paper we consider the evolution of inhomogeneous planetesimal disks caused by in- 
teractions between planetesimals using conventional methods of statistical mechanics (Lifshitz & 
Pitaevskii 1981). The effects of the planetary embryos on the disk evolution will be investigated 
later in Rafikov (2002; hereafter Paper II). We will concentrate on one of the most important pro- 
cesses going on in these disks — mutual gravitational scattering of planetesimals — although our 
rather general approach allows one to treat other relevant phenomena as well. This process can 
proceed in two distinct regimes depending on the amplitude of planetesimal random motion. The 
gravitational attraction between two planetesimals with masses mi and m,2 becomes stronger than 
the tidal field of the central star of mass M c when their mutual separation is less than their Hill 
(or tidal, or Roche) radius rjj, defined as 



where ao is the distance from the central star. When the random velocities of the epicyclic motion 



at ao) their relative approach velocities are small and close interactions can lead to a considerable 
change of the orbital elements of planetesimals. This velocity regime is called shear-dominated 
(or "cold"). It should be contrasted with the other extreme — the so called dispersion-dominated 
("hot") regime which occurs when planetesimal velocity dispersions are bigger than ~ f2r#. In 
this latter case scattering is typically weak which often allows analytical treatment of this velocity 
regime. 

The development of planetesimal disk inhomogeneities driven by a protoplanetary embryo was 
explored in Rafikov (2001; hereafter Paper I) assuming that shear-dominated scattering of planetes- 
imals prevails. It was also assumed in this study that the dynamical properties of planetesimals do 
not evolve as a result of scattering and that the disk always stays dynamically cold. Planetesimal- 




(1) 



of interacting planetesimals are smaller than ~ (^ = \J GM c /clq is the disk orbital frequency 
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planetesimal interactions play the role of effective viscosity in the disk, and tend to homogenize 
it and close up any gap. Nevertheless, this study demonstrated that gap formation is the natural 
outcome of the embryo-planetesimal interaction when the embryo is massive enough. These inter- 
actions were local in character because in the shear-dominated case only planetesimals on orbits 
separated by no more than several r# (corresponding to their encounter) were able to approach 
each other closely. 

It is however more likely that the planetesimal-planetesimal gravitational scattering in realistic 
protoplanetary disks proceeds in the dispersion-dominated (rather than shear-dominated) regime, 
at least on the late stages of disk evolution (see §3). In this case the evolution of planetesimal 
random motion can strongly affect the growth rate of protoplanetary embryos. It is also tightly 
coupled to the evolution of spatial distribution of planetesimals because any change of the energy 
of epicyclic motion comes at the expense of the orbital energy of planetesimals. 

The treatment of the dispersion-dominated case is complicated by the fact that planetesimals 
in this regime can explore different regions of the disk in the course of their epicyclic motion. This 
makes disk evolution a nonlocal process. On the other hand, as we have said before there are natural 
simplifications which are valid in the dispersion-dominated regime. These include the two-body 
scattering approximation (relative velocities are high), Fokker-Planck type expansions (scattering 
is weak), and so on. 

Our paper is organized as follows. In §2 we derive general equations for the planetesimal 
surface density and velocity dispersion evolution. We do this using the Hill approximation which 
is briefly described in §2.1. A Fokker-Planck expansion of the evolution equations, valid in the 
dispersion-dominated regime, is performed in §3. In §4 we derive analytical expressions for the 
scattering coefficients used in these equations. We conclude with a brief summary of our results in 
§5. Some technical details of the calculations and derivations can be found in appendices. 

2. Equations of evolution of planetesimal disk properties. 

Throughout the paper we will assume disk to be axisymmetric. All the following calculations 
also assume a disk with a Keplerian rotation law, although they could be easily extended to the 
case of an arbitrary rotation law. 

As we have mentioned before planetesimal disks must have contained a huge number of con- 
stituent bodies with a wide range of masses. Thus we use a kinetic approach and characterize the 
state of the disk at any moment of time t by a set of three averaged parameters for each mass popu- 
lation: the surface number density of planetesimals and their velocity dispersions in the horizontal 
and vertical directions at every point r in the disk. 

When the orbit's eccentricity and inclination are small, it is convenient to work in the guiding 
center approximation, when the elliptic motion is represented as an epicyclic motion of the body 
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on a small ellipse in horizontal plane with the radial dimension ae (and oscillatory motion in the 
vertical direction with the amplitude ai); the center of this ellipse performs circular motion with 
a semimajor axis a. In the case of Keplerian rotation law frequencies of all these motions are the 
same and equal to the angular frequency of Keplerian motion. In this approximation the roles of 
longitude of ascending node and argument of pericentre are played by two constant phase angles u> 
and r which characterize the position of body on its epicycle at a given moment of time. Random 
motion is defined as the motion relative to the circular orbit of the guiding center, i.e. is represented 
by the epicyclic motion. Eccentricity and inclination thus measure the magnitude of the random 
motion and we will take them to represent planetesimaPs kinematic properties. 

The spatial distribution of planetesimals can be characterized by two kinds of surface den- 
sity. One is the surface density of guiding centers of planetesimals N(a) and is defined such that 
N(a)2irada is the number of planetesimals with guiding centers between a and a + da. It is a 
function of planetesimal semimajor axes only. 

Another parametrization is the instantaneous surface density N mst (r), defined such that 
N mst {r)2irrdr is the number of planetesimals within the disk surface element 2nrdr (we will denote 
instantaneous distance from the central star as r to distinguish it from a). The relation between 
these two surface densities can only be computed if the distribution of planetesimal eccentricities is 
known (we calculate an appropriate conversion between N and N mst in Appendix A). These sur- 
face densities coincide only when planetesimals are on circular orbits (which was the case in Paper 
I). We will mostly work with the surface density of the guiding centers and not the instantaneous 
surface density. 



In the epicyclic approximation, which is valid when e <C 1 and i <C 1, it is convenient to 
introduce the Cartesian coordinate system x, y, z with axes in the r, ip, z directions of the corre- 
sponding cylindrical system and origin at some reference distance ao from the central star. This 
system rotates around the central star with angular velocity Qq = ^(«o)- I n these coordinates the 
unperturbed motion of the body is given by (Goldreich & Lynden-Bell 1965; Henon & Petit 1986) 



where r and to are the constant phase angles we have mentioned before and A is the origin of time 
(it is usually unimportant and can be set equal to 0). Dimensionless relative semimajor axis a of 
the body (the location of its guiding center) h is defined as 



2.1. Hill equations. 



x = hao — eao cos(ftot + r) + 0(e 2 ,i 2 ), 
3 

y = —-hao(Qot — A) + 2eao sin(f2o£ + t) + 0(e 2 ,i 2 
z = iao sin(f2ot + oj) + 0(e 2 , i 2 ), 



) 



(2) 



(3) 
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Instead of the pairs e, r and i, uj it is sometimes convenient to use the eccentricity and inclination 
vectors e and i: 

e = (e x , e y ) = (ecosr, esinr), 

i = (ix,iy) = (icosu, ismu). (4) 

We can also use other natural simplifications. One of them is motivated by the assumption 
that the masses of planetesimals rrij are much smaller than the mass of the central star M c — a 
condition which is always fulfilled in planetesimal disks. In this case it was demonstrated by Henon 
& Petit (1986) that the equations describing the 3-body interaction (star and planetesimals 1 and 
2) can be significantly simplified. In particular they have demonstrated that relative motion of 
interacting bodies (r = ri — T2) and the motion of their barycenter [r^ = (?ni r i + wi2 r 2) / + ^2)] 
separate from each other. In the course of 3-body interactions the barycenter motion is not affected 
at all and remains in unperturbed epicyclic motion at all times. 

We define the relative Hill coordinates of the two interacting bodies (x,y,z) with masses 
mi *C M c and mi <C M c as 

x = (xi - x 2 )/r H , y = (yi- yi)/rH, z = (zi- z 2 )/r H . (5) 

Here r# is the Hill radius defined by equation (1). Instead of m we will often use n = m/M c [with 
this notation the Hill radius can be written as rn = oo(a*i + A*2) 1 ^ 3 ]- Note that the Hill radius is 
only meaningful for a pair of bodies and its definition is similar to the definition of tidal radius. We 
will also define reduced eccentricity and inclination vectors e = ao(ei — e2)/r#,i = oo(ii — \2)/th 
and the reduced impact parameter h = ao(hi — h2)/rn- 

With these definitions the equations of relative motion can be written as (Henon & Petit 1986; 
Hasegawa k, Nakazawa 1990) 

d 2 x _ dy_ _ _ d(j) p 
dT 2 dT 6X ~ dx 
d 2 y dx_ ^ (typ 
dT 2 dT~ dy' 
d 2 z d(f> p 
dT 2 + z= ~qz" 

where 4> p is a (dimensionless) potential of interaction between the two bodies and T = flot. For 
the Newtonian gravitational potential between planetesimals 

cj )p ( Y )=-(x 2 + y 2 + z 2 )- 1 / 2 . (9) 

An important feature of these equations is that they do not contain any parameters of the interacting 
bodies (such as their masses, etc.). All the physically important information is embedded into the 
definitions of Hill radius and relative coordinates. Also, the outcome of the interaction between 



(6) 
(7) 
(8) 
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the two bodies depends only on their relative orbital parameters. This universality is a very useful 
property which we will widely exploit later on. Equations (6)-(8) possess an integral of motion - 
the Jacobi constant 



-l§) 2 +(f) 2 + (§) 2 -3^ + %. (10 

Conservation of this quantity poses an important constraint on the changes of the orbital elements 
of interacting bodies in the course of their gravitational interaction. 

One can easily see that when the interaction is absent equations (2) represent a solution of the 
system (6)- (8) if we replace all quantities in (2) by their relative (and scaled by r#) analogs. But 
in the presence of the interaction, orbital parameters are no longer constants and will evolve with 
time. However, we can still represent epicyclic motion by equations (2) keeping in mind that e, i, h, 
etc. are now instantaneous osculating values of the orbital elements. Their evolution is governed 
by the following set of equations (Goldreich & Tremaine 1980; Hasegawa & Nakazawa 1990): 

dh _ d(f>p dX _ 4 d(f) p 2 .d<f>„ 
- ~ 2 ^?r> - ttt^t" - ri J ~ A )- 



dT dy ' dT 3h dx h" ' dy ' 

de x . d<t> p d(j) p de y d(p p . d</> p 

It = " smT ^ " 2cosT W ~dr = " cosT ^ + 2smT W 

§ = -cos7* § = sinT^. (11) 

dT dz dT dz y ' 

The relative velocity r of planetesimal motion in Hill coordinates (dimensional) is given by [see 

(2)] 

3 ~ 

x = eru sin(T + r), y = — i^ ir H + 2er# cos(T + r), z = iru cos(T + to). (12) 

In many applications it is better to use the velocity defined relative to the local velocity of the 
circular motion. Clearly this velocity v is related to f as 

3 

v = r+-xn y , (13) 
where n y is a unit vector in the y-direction. This definition implies that 

v x = x = eru sin(T + r), v y = y + -x = cos(T + r), v z = z = iru cos(T + u). (14) 

With the use of equations (2) and (14) one can express e, i and h as functions of r and v which will 
later be used in §4. Using (2) and (12) we can rewrite the expression (10) for the Jacobi constant 
in the following form: 

J = e 2 + i 2 - H 2 + 2(f) p (15) 

(our definition of the Jacobi constant may differ from others by a constant additive and/or multi- 
plicative factors, cf. Hasegawa & Nakazawa 1990). 
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2.2. Derivation of the master equation for the distribution function. 

Based on the previous discussion it is natural to describe the position and velocity state of 
a planetesimal with mass rrij by its guiding center radius hj (relative to the reference radius an), 
vector eccentricity ej and inclination ij 1 (note that this is not a canonical set of variables!). The five- 
dimensional vector Tj = (hj,ej, ij) will be used to concisely denote this state and dTj = dhjdejdij 
will denote a phase space volume element around Tj. 

We set N(m, h, t)dm to be the dimensionless surface number density of planetesimal guiding 
centers in the mass interval (m,m + dm). It differs from its dimensional analog N(a) by a factor 
of <2q. We introduce the notion of a planetesimal distribution function (P1DF) f(m,T,t) such that 
2irf(m,T,t)dTdm is the number of planetesimals in phase space volume dT and in mass interval 
(m,m + dm). Clearly, 

J J f(m,T,t)dedi = N(m,h,t), (16) 

where the integration is performed over the 4-dimensional space of vector eccentricity. For the 
purpose of brevity we will further use fj(T,t) and Nj(h,t) instead of f(m,j,Tj,t) and N(rrij,h,t). 

In what follows we will consider the distribution function of planetesimal vector eccentricities 
and inclinations to be Gaussian (or Rayleigh, or Schwarzschild). In other words, we assume that 
fj(T)dT = N(h)ipj(e, i, h)dT, where 



ipj(e, i, h) = exp 



e 2 i 2 



2a%{h) 2a^(h) 



(17) 



Here a e j,aij are the r.m.s. eccentricities and inclinations of planetesimals of mass rtij at point h, 
given by 

a ej = (e 2 ) 1 / 2 /^, a l3 = (i]) 1 ' 2 /^. (18) 

Greenzweig Sz Lissauer (1992) and Ida & Makino (1992) have demonstrated that this assumption of 
Gaussian distribution is very good in their N-body simulations for large values of planetesimal ec- 
centricities and inclinations [see also Ohtsuki (1999) for the discussion of the validity of distribution 
(17) in circumplanetary disks]. 

Consider now planetesimals with masses mi (type 1) and m,2 (type 2) (we usually denote the 
particle under consideration by subscript 1 and scattered particles by subscript 2). Scattering can 
bring planetesimals of type 1 into the phase volume element dT and can also scatter particles out 
of it. The number of encounters in time dt between particles of type 1 initially in the phase volume 
element dT a and those of type 2 initially in dTf, is 

dN coll = u(r a , r 6 )/i(r )/ 2 (r 6 )dr dr 6 di, (19) 



1 These variables are dimensionless but not scaled by m as in §2.1. 
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where v(T a ,T b ) is the rate of encounters. 

Let us introduce the differential probability P(T a ,Tb, AT a ) such that P(T a ,Tb, AT a )d(AT a ) 
is the probability that particle 1 changes its state from T a to the region dAT a around T a + Ar a in 
an encounter with particle of type 2, if before the collision they were in states r a , T b respectively. 
We can then write down the evolution equation for P1DF of particles of type 1: 

oo 

dfi(T) +t dh(T) 



dt dT 

o 



= J dm 2 J dT a dT b 



V (r a , r 6 )/i(r a )/ 2 (r 6 )P(r , r fe , r - r a ) - i/(r, r 6 )/i(r)/ 2 (r 6 )P(r, r 6) r a - r) 



(20) 



Here the second term in the l.h.s. represents evolution caused by reasons other than mutual 
gravitational scattering (e.g. gas drag, migration, non-Keplerian gravitational forces, etc.). The 
r.h.s. of (20) is the collision integral for scattering between planetesimals; its first term represents 
particles entering dT, while the second one represents planetesimals leaving dr. This is the most 
general form of evolution equation, which will be further simplified using additional assumptions. 



2.3. Local approximation. 



When the masses of the interacting bodies are much smaller than M c their mutual Hill radius 
is much smaller than the distance to the central star do- We will also assume that the sizes of 
epicycles of individual bodies are small compared to ao- As a result, we can safely use the Hill 
formalism (§2.1) to describe planetesimal scattering. 

The local character of the interaction means that the encounter rate u(T a ,T b ) is determined 
by the local shear and can be expressed as 



u(T a ,T b ) = Mh a ) - n(h b )\ = 2\A\\h a - h b \, 



(21) 



where A = (r/2)(dfi/dr) is a function characterizing shear in the disk (Oort's A constant), A 
— (3/4)f2o for a Keplerian rotation law. Using this expansion we can rewrite equation (20) as 



dh(T) +t dh(T) 



dt 



dT 



2\A\ J dm 2 J dT a dT b 



\h a - h b \h{T a )h{T b )P{T a ,T b ,T - T a ) - \h - / lfc |/ 1 (r)/ 2 (r fe )p(r, r b , r a 

We can easily integrate the second term over dT a because it is clear that 

J dT a P(T, T b , r a — r) = 1, 



(22) 



(23) 
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since this integral represents the probability of scattering planetesimal 1 (initially at T) anywhere 
in the course of interaction with planetesimal 2 (initially at T b ). Then we can also integrate f2{^2) 
over de2di2 to find 

/roo 
dT a dT b \h- h b \f 1 (T)f 2 (T b )P(T,T b ,T a -T) = / dh b \h - h b \h{T)N 2 {h b ). (24) 
•J — oo 

To deal with the first term in the r.h.s. of (22) we first recall that in the Hill approximation 
all scattering properties depend only on relative coordinates and velocities of planetesimals. This 
means that 

p(r a ,r 6 ,Ar a ) = p(r a -r 6 ,Ar a ), (25) 

Now, let us introduce the new relative variable 

f r = (h, e r , i r ) = - — — ^ = ■ — ^ (h a - h b , e a - e b , i a - i b ) , (26) 

where T a and r a characterize the orbital elements of particles 1 and 2 correspondingly The change 
of T r in the course of an encounter is AT r , and depends only on r a and T b through the combination 
r r . Then, bearing in mind conservation of the center of mass coordinate miT a + iri2T b , we obtain 
that r a , T b can be expressed in terms of T [in the first term of (22) this is the value of r a after the 
encounter], T r , and AT r as 

r a = r - — ^—(jii + ^ 2 ) 1/3 Afv, r b = r - ( Ml + ^) 1/3 ( — ^— Af r + r r ) . (27) 

Instead of the probability P(T a , T b , V — T a )d (AT a ) we will work with the probability of scat- 
tering in relative coordinates P r (T r , AT r )d(AT r ). From the equation (27) one can easily find that 

P r (f r , At r ) = ^ P(T a - T b , AT a ). (28) 

(m + fi 2 ) w/6 

This new probability function is independent of the masses of the interacting particles. It also 
follows from the properties of the Hill equations that 

P r (f r , Af r ) = 6 [AT r - Af ac (f r )] , (29) 

i.e. the initial state of the scattering planetesimals unambiguously determines their final state 
through the scattering function Ar sc (r r ). 

Using these definitions we can finally put the P1DF evolution equation (22) in the following 
form: 

dh(T) , dTdh(T) 



+ f ^§r^ = _2|A| / dm2{f11 + ^ 2)2/3 / m 



dt 

o 



/i(r)iv 2 (/ i -(/x 1 + // 2 ) 1 /3^) (30) 

(^ + /x 2 ) 4/3 y de r d\ r d(At r )h(T a )f 2 (T b )P r (t r ,AT r 
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with r, r a ,rfc,rY and Ar r related by equations (26) and (27). 

We will use this general equation for the evolution of the P1DF to obtain formulae describing 
the behavior of the spatial and kinematic properties of planetesimal disk. 



2.4. Evolution of surface density. 



In this paper we will be mostly interested in the evolution of the planetesimal disk due to 
the gravitational scattering between planetesimals. Thus, we set the second term in the l.h.s. of 
equation (30) to (it can be easily reinstated when the need arises). At this stage we also replace 
subscripts "a" and "b" with "1" and "2" since now r a and have the only meaning of initial 
orbital elements of particles 1 and 2. 

To obtain the equation of evolution of Ni(h,t) let us multiply equation (30) by dedi and 
integrate. Using equation (16) we can easily do this with the l.h.s. and the first term of the r.h.s. 
To treat the second term in the r.h.s. of (30) we substitute the distribution function (17) into 
equation (30) and integrate over de using (27) and the following identity: 



/ 



de 

~a 9 2 2~~ eX P 



2**2 J 



2vr (a 2 el + a 2 e2 ) 



cxp 



2 + <? 2 e2 ) 



ip r (e r ) , (31) 



with ip r (e r ) being the distribution function of relative eccentricity as defined by (31). An analogous 
identity defines the distribution function of relative inclination ij) r (i r ). These identities show that 
the distribution function for the relative motion of planetesimals with r.m.s. eccentricities and 
inclinations a e i,an and cr e2 ,ai 2 is also Gaussian, with dispersions of the relative eccentricities and 
inclinations given by 

(32) 



^-2 , „2 



^ 2 , „2 

a n + ca- 



using this result we are able to write the following equation for the evolution of N\{h): 

OO OO f 

0N\ (h) . .. r r ~.~ I 



where 



dt 



= -2\A\ J dm 2 {m + fi 2 ) 2/3 J dh\h\lN 1 (h)N 2 [h-(fi 1 + fi 2 ) 1 / 3 h] 



- (m + V2) 4/3 J de r di r d(Af r )V'r(e r ,a er )V'r(ir,^ r )^i(/ i i)A^ 2 (/i 2 )P r (f r ,Af r )| , (33) 



<? 2 er {hi,h 2 ) = a 2 el {hi) +<J 2 e2 {h 2 ), a 2 r {hi,h 2 ) = a^hi) +a 2 2 (h 2 ), 
h 1 = h + D(Ah), h 2 = h + D(Ah) - ( Pl + fi 2 ) 1/3 h, 
with D(Ah) = — ^——Ah. 



(/il+/U 2 ) 2 /3 



(34) 
(35) 
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In the last term of (33) one can perform the integration over d(AT r ) using equation (29) [which 
reduces to simple replacement of Ah by Ah sc (h, e r , i r ), see equation (29)]. Also it is more natural 
to use the reduced velocity distribution function ip r (e r , i r ) defined as 



Vv(e r ,i r ) = (m + / u 2 ) 4/3 Vv(e r )Vv(ir) 



Air 2 a 2 a 2 



cxp 



el 



where 



0V. 



e * (M1+M2) 1 / 3 ' 

Using all these results and definitions one can finally transform equation (33) into 



aZVi(h) 
dt 



-2\A\ j dm 2 ( / u 1 + /t / 2 ) 2 / 3 J dh\h\\ N^hWh- (pi + ^ 2 ) 1/3 h] 

-oo ^ 

- J de r di r 4> r (e r X)Ni[h 4 D(Ah sc )]N 2 [h 4 D(Ah sc ) - (m 4 /x 2 ) 1/3 /t]| . 



(36) 



(37) 



(38) 



The functions iV 12 (/i) and ip r (e r ,i r ) are also functions of time. This equation is analogous to the 
evolution equation (6) of Paper I but now it can also take into account the dependence of scattering 
properties on the random motion of interacting bodies. 



2.5. Evolution of random velocities. 

Let us multiply equation (30) by e 2 <iedi and, again, integrate over the whole vector eccentricity 
and inclination space. Since integrating the l.h.s. of (30) and the first term in its r.h.s. is again 
trivial [by definition J e 2 fj(e, i)dedi = 2Nja 2 ], we concentrate on the second term in the r.h.s. We 
notice using (31) that the integration of ipi^e 2 de reduces to 



cxp 



2 (<& + <£) 



e 2 de 



4vrV 2 - 2 



■ cxp 



= Vv (e r ) 



2 4 



rr 4 p 2 
a el e r 



4 



,:2J 



°ll + ^e2 
2 



M2 



A*i + M2 



-Ae r 



cr el e r 



M2 



Ml + M2 



(Ae r ) 2 4 2 



M2 



0"; 



el 



Ml + M2 <i + cr, 



e2 



2 * 

e2 



(39) 



[Integration over cZi yields ip r (i r ), see equations (31)]. Using again the P1DF in the form (29) we 
find that integration of the second part of the r.h.s. of (30) gives 



J de r di r 4> r (e r , \ r )Ni[h + D(Ah sc )]N 2 [h + D(Ah sc ) - (jn + fi 2 ) 1/3 h] 



a 2 ! a 2 n 

2 el ez 



rr A p 2 



4 



*fl + ^e2 (^el + CJ e 2 ) 2 VMl + M2 



M2 



(Ae sc ) 2 4 2 



M2 



cr 



el 



Ml + M2 cr^ 4 rr e2 



-e r • Ae s 



,(40) 
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where we have again used function ip r (e r ,i r ) defined in equation (36), and 

Ae sc (h, e r , i r ) = (m + ^ 2 ) 1/3 Ae sc (h, e r , i r ) 
comes from equation (29). 



(41) 



Using these results we can rewrite the evolution equation of the r.m.s. eccentricity of planetes- 
imals of mass m\ in the following form: 



d_ 

at 



-2\A\ J dm 2 (fii + n 2 ) 2/3 
o 



dh\h\ 



xi 2a 2 el {h)N 1 {h)N 2 [h - (m + fi 2 ) 1/3 h] - J de r d\ r ^ r {e r , I T -)iVi(^i)7V 2 (/i 2 ) 



+ 



n A P 2 

a el e r 



+ 



/*2 



(Ae sc ) 2 + 2 



M2 



+ (P'el + fT e 2 ) 2 \Ml + ^2 

As usual, it is implied in the last term of this equation that 



Ml + M2 cj^ + cr 2 2 



Ae s 



(42) 



°li = a ei i h i) i °"e2 = °"ei i h 2) , where now 

/ii = + D(Ah sc ), h 2 = h + D{Ah sc ) - (m + /i 2 ) 1/3 /i, 



(43) 



i.e. only the kinematic characteristics of planetesimals of different masses at their locations before 
the encounter are important. 

It should be remembered that the integral over de r di r in the r.h.s. of equation (42) cannot 
be easily evaluated in general, because Ah sc and Ae sc depend on both e r and i r . However, if 
the disk kinematic properties and surface density are homogeneous (i.e. a e j, a^j, and Nj are 
independent of h), then the integration over de r di r can be easily performed on the first two terms 
in the square brackets of equation (42) and they vanish. Planetesimal eccentricities would then 
evolve only because of the presence of the last two terms in square brackets in (42). This leads 
us to the conclusion that these nonvanishing terms [third and fourth in the r.h.s. of (42)] are 
responsible for the planetesimal self-stirring — kinematic heating in the course of gravitational 
scattering of planetesimals (which is present even in a completely homogeneous disk) . The equation 
of eccentricity evolution in homogeneous disk then assumes the form [using (42)] 



dt 



oo 

= \A\ J dm 2 {^ + ^ 2 ) 2 ^N 2 



M2 



Mi + M2 



ill + 2 



M2 



el 



Ml + M2 cr 2 el + a 2 e2 



Ho 



(44) 



where self-stirring coefficients are defined as 

oo oo 

Hi = j dh\h\ j de r di r -i/v(e r , i r ) (Ae sc ) 2 , H 2 = j dh\h\ j de r di r tj) r (e r , i r )(e r • Ae sc ). (45) 
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Equation (44) and definitions (45) are analogous to equations of evolution of planetesimal eccen- 
tricity dispersion derived by other authors for the case of a homogeneous planetesimal disk (Ida 
1990; Wetherill & Stewart 1993; SI). 

The other terms in the r.h.s. of equation (42), which disappear in the homogeneous disk, 
must describe the transport of horizontal energy due to gravitational scattering, and are caused by 
(1) the gradients of r.m.s. eccentricity and inclination of planetesimals, and (2) the gradients of 
planetesimal number density To the best of our knowledge, this part of eccentricity evolution has 
never been taken explicitly into account before. In the general case of an inhomogeneous disk (e.g. 
in the presence of a massive embryo which can naturally induce gradients of disk properties by its 
gravity) both contributions (self-heating and transport) are important. 

In analogous fashion we can write the equation for the inclination evolution: 

oo oo 

^ [ZN^alih)} = -2\A\ f dm 2 (fi 1 + f i 2 ) 2 ^ J dh\h\ 

-oo 
{ 2al(h)N 1 (h)N 2 [h - ( Ml + fl 2 ) 1 / 3 h] - J dMir^r(^,ir)Wi(fci)AT 2 (/l 2 ) 

2 44 . °m j ^2 y (A . )2 2 _^ 0%_ 

4+4 (4+4) \Ml+M2/ /ii + /i2 Ofi + <7f 2 



(46) 



[the definition of Ai sc is analogous to definition of Ae sc in equation (41)]. We cannot simplify 
these equations further without making additional assumptions about the scattering properties of 
planetesimals. These additional approximations will be made in the next section. 



3. Fokker-Planck expansion. 

Significant simplification can be achieved if gravitational scattering is weak, i.e. changes in 
quantities characterizing the planetesimal state T are small compared with the average values of 
these quantities before the encounter. Then a Fokker-Planck type expansion can be performed 
(Lifshitz k Pitaevskii 1981; Binney & Tremaine 1987). 

This situation is usually the case for planetesimal-planetesimal interactions because the evo- 
lution of the disk's kinematic properties prior to the emergence of embryos would likely lead to a 
considerable dynamical "heating" of planetesimal population. This heating brings planetesimal- 
planetesimal interactions into the dispersion-dominated regime when two interacting planetesimals 
of types 1 and 2 have a\ r + af r » 1. Then encounters occur at high relative velocities and large- 
angle scattering is rare. For example, a rocky planetesimal with the radius 50 km at 1 AU has 
a corresponding mass of ~ 10 21 g and Hill radius ~ 10~ 4 AU. The critical velocity determining 
the lower boundary of the dispersion-dominated regime for the gravitational interactions between 
such planetesimals is about 3 m s" 1 , which is likely to be exceeded by the velocity of planetesimal 
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epicyclic motion (e.g. see the results of coagulation simulations by Inaba et al. 2001). Thus, the 
assumption that planetesimal-planetesimal scattering occurs in the dispersion-dominated regime is 
usually justified in protoplanetary disks. 

The necessary condition for the Fokker-Planck expansion is that changes of h,e r ,i r are small 
compared with their typical magnitudes. For example, for the dimensionless impact parameter h 
this means that aAh should be much smaller than the typical length scale on which the surface 
density of planetesimals varies. However, this approximation does not require that the relative 
distance between interacting bodies ah must be small in comparison to this typical scale. This is 
in contrast to the local approximation that is often used along with the Fokker-Planck approach. 
Indeed, when the planetesimal disk is dynamically hot, every planetesimal can interact with all 
others within the reach of its own epicyclic excursion. Since eccentricity is large in this velocity 
regime the semimajor axis separation of interacting planetesimals ~ ae could be of the order of or 
bigger than m\ then local approximation (which would assume that semimajor axis separation of 
interacting planetesimals is small) becomes inapplicable. This complication is explicitly taken into 
account in our further consideration. 



3.1. Surface density evolution equation. 



To implement the Fokker-Planck treatment of the surface density evolution equation (38), we 
first expand the second term of the r.h.s. of equation (38) in a Taylor series in D(Ah sc ) up to the 
second order. Using the notation a = (/ii + [i 2 ) 1 ^ and h — (hi + jx^^h = h — ah we obtain that: 



/ 



de r d\ r ^ r {e r X)Ni[h + D(Ah sc )]N 2 [h + D(Ah sc ) - (m + fi 2 ) 1/3 h] 

IJ-2 d 



+ 



N 1 (h)N 2 (h - ah) 

1 y 2 d 2 

2 (/xi + /x 2 ) 4 / 3 dh 2 L 



(/ii+^2) 2 / 3 dh 
{(Ah) 2 )N 1 (h)N 2 (h - ah) 



{Ah)Ni{h)N 2 {h- ah) 



(47) 



In obtaining this expression we have taken into account that the distribution function ?/v is a 
function of a 2 r , af r , and, thus, should also be expanded in D(Ah sc ). In equation (47) the moments 
of Ah are 



{(Ah) 



I 



Ah sc (h,e r ,i r ) ip r (e r ,i r ) (48) 
where = 1,2, and ((Ahf) is a function of h,cr er ,Oi r . Substituting (47) into (38) we obtain 

dh 2 

where 



(49) 



00 00 
(h)=2\A\ j dm 2 fi 2 J dh\h\(A~h)N 2 (h- ah), 



(50) 
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T?(h) = \A\ [ dm 2 - [ dh\h\((Ah) 2 )N 2 (h-ah). (51) 

-oo 

Thus, instead of the integro-differential equation (38) we have obtained the partial differential 
equation (49). 

The integration over dh in definitions of transport coefficients and takes into account 
the non-zero range of the planetesimal disk over which a given planetesimal can experience encoun- 
ters with other bodies (nonlocal scattering). In the local approximation, (Ah) and ((Ah) 2 } would 
fall off rapidly with increasing h; expanding the integrand in the expression for in h one could 
easily derive the local Fokker-Planck evolution equation. This equation has already been discussed 
by Petit & Henon (1987b) and in Paper I. 



3.2. Random velocity evolution equation. 

Now we perform the Fokker-Planck expansion of equation (42) for the velocity dispersion. 
Again, we need to concentrate on the expansion of the terms under the de r di r integral with respect 
to D(Ah sc ). 

The Fokker-Planck expansion usually assumes keeping only the terms of the two leading orders: 
the first one is linear in the expansion parameter, but can suffer strong cancellation effects when 
averaging over the scattering outcomes is performed. The second term is quadratic but it does 
not suffer from cancellation during the averaging and in general has magnitude similar to the first 
one. In our case terms such as (Ah sc ) 2 e ■ Ae sc or Ah sc (Ae sc ) 2 are third order in the perturbation 
and should be neglected (we will comment on the order of their smallness in §4). But terms like 
Ah sc e ■ Ae sc or (Ah sc ) 2 e 2 should be kept in our expansion. 

Performing this procedure as in §3.1, after some cumbersome but straightforward transforma- 
tions, we can write the Fokker-Planck equation for eccentricity evolution in the following form: 

| [2^] = r* Nl - J- (TfiVi) + (T|i\q) (52) 



where 



T e (h) = 2\A\ J dm 2 (m + ^ 2 ) 4 / 3 J dh\h\N 2 (h - ah) 

-oo 

-^-V <(Ae) 2 > + 2 ^ 2 ./ gl . 2 (e • Ae> 

oo oo 

Tf (h) = 2\A\ J dm 2 /x 2 (/xi + n 2 ) 2/3 J dh\~h\N 2 (h - ah) 



(53) 
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2 ° eia % (Ah) + — 2 a * (e 2 Ah) + 2 m _ 2 ael _ 2 ((§ • Ae)A/i> 

OO OO 

T e 2 (h) = \A\ J dm 2t ? 2 j dh\h\N 2 (h - ah) 



(54) 



x 



2#^| r <(A/ i ) 2 > + gLg (e 2 (Afr) 2 ) 



(55) 



The new scattering coefficients used in these equations are 

((Ae) 2 ) = y 7p r (e r , i r ) (Ae sc ) 2 <ie r <ii r , (e • Ae) = ^ ^ r (e r , i r )e r • Ae sc <ie r di r , 
(e 2 Ah) = J ^ r (e r ,i r )e 2 ,Ah sc de r di r , (e 2 (Ah) 2 ) = J tj) r (e r ,i r )e 2 (Ah sc ) 2 de r di r , (56) 
((e • Ae)Ah) = J Vv(e r , i r )A^, sc (e r • Ae sc )de r di r , 

These coefficients, like the old ones ((Ah), ((Ah) 2 )) are functions of h, a^, a 2 2 , a 2 l} a 2 2 . Analytical 
expressions for these coefficients in the two-body approximation will be derived in §4. An analogous 
evolution equation can be written for a\. 

The first term in the r.h.s. of equation (52) is responsible for the gravitational stirring of 
eccentricity, while the last two terms describe energy transport and disappear in homogeneous 
planetesimal disks. The coefficients in (49) and (52) are nonlocal quantities and this is reflected in 
the presence of an integration over dh in their definitions. 

Equations (49)-(51) and (52)-(55) describe self-consistently the evolution of the surface den- 
sity, eccentricity, and inclination at every point of an inhomogeneous planetesimal disk due to 
planetesimal-planetesimal gravitational encounters in the dispersion-dominated regime. The prin- 
cipal approximation is that we follow moments rather than whole distribution function. To close 
this system we only need to supplement it with the expressions for various scattering coefficients 
appearing in evolution equations. 

4. Scattering coefficients in the dispersion-dominated regime. 

In this section we derive analytical expressions for the scattering coefficients used in equations 
(49)-(51) and (52)-(55). In so doing we will always assume that interaction between planetesimals 
is in the dispersion-dominated regime, for the reasons given at the beginning of §3. In the case of 
embryo-planetesimal scattering it is less clear that the scattering is dispersion-dominated since the 
embryo mass and Hill radius are much larger; we consider these additional details in Paper II. 

In the dispersion-dominated regime most interactions occur between planetesimals separated 
by \h\ » 1 (or \h\ » th)- In this case they can be divided into those that can experience close 
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Fig. 1. — Separation of different scattering regimes in the e — i phase space. 

encounters, and those which do not approach each other closely. For the planetesimals of the 
first type most of the change of their orbital elements occurs near the closest approach point. 
Moreover, this change occurs so rapidly that one can consider scattering to be instantaneous and 
local, i.e. to occur at some point and not over an extended part of the trajectory. These observations 
considerably simplify the analytical treatment of the problem because one only needs to consider a 
small region near the encounter point. Clearly, only planetesimals with e > \h\ fall in this category 
(all the orbital parameters are relative and we drop subscript "r" in this section). It will turn 
out that planetesimals in this part of the parameter space (e > \h\ and arbitrary i) are the most 
important contributors to the scattering coefficients. We will call this part of the e— i space "Region 
1" and devote §4.1 to its study (see Figure 1). 

Some planetesimals cannot approach closely — those with e < \h\. This part of the phase 
space can be split into two more regions: one with i < \h\, we will call it "Region 2" (see Figure 
1), and the other with i > \h\ ("Region 3"). In the Region 2 planetesimals accumulate changes of 
their orbital parameters along a significant portion of their orbits, which stretches approximately 
by several times h along the y-direction. In a zeroth approximation, the planetesimal disk in this 
part of phase space can be treated as not possessing random motion at all. Encounters between 
planetesimals in Region 2 lead to smaller change in the orbital elements than encounters in Region 
1 because close encounters cannot take place in Region 2 (Hasegawa & Nakazawa 1990; SI). In 
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Region 3 typical distances at which interaction occurs are similar to those of Region 2 but relative 
velocities are higher (they are magnified by the large value of i), which makes scattering in Region 
3 weaker than either Region 1 or Region 2. 

In general all 3 regions will contribute to the scattering coefficients for a given value of h. 
However, their contributions are not of the same order and their relative magnitudes strongly 
depend on the value of h. Assume first that \h\ < a e . Then the distribution (17) implies that 
most of the planetesimals belong to Region 1 (where e > \h\). Since every individual scattering 
event in Region 1 is also stronger than those happening in Regions 2 and 3 we can conclude that 
Region 1 strongly dominates the scattering coefficients when \h\ < a e . However, if a e < \h\ the 
number of planetesimals in Region 1 of the phase space becomes exponentially small, so that distant 
encounters (Region 2) will dominate the scattering coefficients. 

For planetesimal-planetesimal scattering we are mostly interested in quantities integrated over 
some part of planetesimal disk, such as the coefficients Tj^ and T e - in equations (49) and (52). They 
are rather insensitive to the details of the distribution of the integrands with h. Moreover, one can 
show that the contribution of distant encounters (Regions 2 and 3) to the integrated quantities is 
still subdominant compared to those from Region 1 (SI). For this reason we will derive only the 
contribution to the scattering coefficients from Region 1 of the phase space and will completely 
neglect distant encounters when discussing the planetesimal-planetesimal interactions. The role of 
distant encounters in the context of the embryo-planetesimal scattering is more important and will 
be discussed in Paper II. 

4.1. Scattering coefficients due to close encounters. 

Since our disk evolution equations use h, e, and i as coordinates and since we are interested 
in dealing with inhomogeneous planetesimal disks it is natural to use the Hill formalism (Ida et 
al. 1993; Tanaka & Ida 1996; SI) to calculate scattering coefficients. This is opposed to the local 
velocity formalism of Hornung et al. (1985) which would be rather awkward for our purposes. 

The treatment of the close encounters in the dispersion-dominated regime is significantly sim- 
plified by the fact that the interaction has a 3-dimensional character — planetesimals can approach 
the scatterer from all directions. This is in contrast to distant encounters which are intrinsically 
2-dimensional. A characteristic feature of 3-dimensional scattering is that if we order planetesi- 
mals by their impact parameter, then equal logarithmic intervals in impact parameter contribute 
equally to the scattering (Binney & Tremaine 1987). In the dispersion-dominated regime most of 
this interval corresponds to minimum separations less than rjj- In this case close encounters can be 
treated by using the two-body approximation which neglects the influence of the third body (Sun) 
on the process of instantaneous local scattering. 

We will derive analytical expressions for the scattering coefficients (48) and (56) averaged over 
dedi. Switching from e and i back to e~,i,T, and to (absolute values of relative eccentricity and 
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inclination and corresponding constant phases, see §2.1) we find that dedi = ede idi dr du. The 
phases t/2tt and uj/2h are distributed uniformly in the interval (0, 1), while the distribution of e, i 
has a form [equation (17)] 



ede idi 
ip[e, i)dedi = ~ 2 ~2 ex P 



<*t<n 



2a 2 2a* 



(57) 



We will perform the averaging in two steps: first we average over phases r and oj (these 
averages will be denoted by (...) Tll) ), and then we average over the absolute values of eccentricity 
and inclination (these final averages will be denoted by simply (...)). 

Thus, initially we fix h, e, and i. The idea of the calculation is to assume the flux and velocity 
of incoming particles near the location of the scatterer (the planetesimal on which we center the 
reference frame of relative motion) to be represented by their unperturbed values and use them as 
an input for the two-body scattering problem. These quantities will be given by their values at the 
scatterer's location r = which means that we only retain zeroth order terms of their expansion 
near r = 0. Then the components of relative planetesimal approach velocity v in the r, (p and z 
directions are given by evaluating the time derivatives of equations (2) at x = y = z = 0: 



= ±Qr H \J e 2 - h 2 , v y = Qr H ^, v z = ±Q.r H i. (58) 

Everywhere in this section velocities v x ,v y ,v z are assumed to be defined relative to the local circular 
orbit (i.e. v y = y — 2 Ax, see §2.1). 

Next we replace the variables r and uj by the impact parameter, I, and the angle 4> of the 
orbital plane (see Appendix 8. A of Binney & Tremaine 1987 2 ). We will use I and 4> as an equivalent 
set of variables and replace averaging over t/2it and lo/2tt by the averaging over I and <fi. Since we 
assume that the planetesimal flux is locally homogeneous around the scatterer the distribution of 
planetesimal trajectories in (p is uniform. 

Let us set To and ojq to be the values of r and uj for which the planetesimal orbit passes through 
the location of the scatterer r = 0. There are 4 such sets of r and uj [because v x and v z have a 
sign ambiguity in (58)]. It was demonstrated by Ida et al. (1993) that the planetesimal trajectories 
have impact parameters smaller than / only when the phases r and uj of their orbits lie within 4 
ellipsoids near 4 sets of To, ujq in r — uj space. The total surface area covered by these ellipsoids is 
given by 

A(l)-- — — 1 

~ 3r 2 H ^r H i\h\^/e 2 -~h 2 

where 



e 2 + i 2 - (3/4) h 2 (60) 



2 In the notation of Binney & Tremaine (1987) I = b. 
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is the magnitude of the planetesimal approach velocity [see (58)]. Thus, particles having impact 
parameters in the range from I to I + Al cover a surface area {dA/dl)Al inr-w phase space. This 
allows us to replace the average of some quantity / with respect to drduj / '(47r 2 ) by an integration 
with respect to dld(j) using the following conversion: 

2-7T 2n Imax 27T 



The upper cutoff distance l ma x will be specified later. To be precise, the transition represented 
by (61) is legitimate only when I <C th- Some planetesimals pass the scatterer at distances > rn 
but they contribute only weakly to the overall change of the orbital elements (because of the 
aforementioned dependence on logarithmic and not linear intervals in I). For this reason the errors 
that arise from using (59) at |r — tq\ ~ 1 or \u — ojo| ~ 1 will hardly affect the integrals in (61). 

Using equations (2), (4), and (14) one finds that the various combinations of the changes of 
the relative orbital parameters that we need can be written as 

A(e 2 ) = [A(vD + 4A(t#] , A(i 2 ) = ^A^ 2 ), (62) 

(the collision is assumed to be instantaneous here). Obviously, A(vf) = 2v{Avi + (Avi) 2 . From 
equations (4) and (14) we can also derive that h = (2v y + fix) / (fir h), meaning that 



Ah = 2^ and (Ah) 2 =4^^. (64) 

We also have 



(e-Ae)A^ = 2 ^ A ^ 3 + r 3 4 ^ (A ^ )2 , (i • Ai)A~h = 2 ^f "» . (65) 

To average these expressions over r, u we first integrate them over d(f>/2ir which means that we 
need to know (vi)^ and (v-iVj)^ (here i = x,y,z and {...)</, means averaging over d(j)/2ir). Two-body 
scattering in terms of coordinates I, <p was considered in detail elsewhere (Binney & Tremaine 1987) 
and we will simply use known results for these averages here: 

V i //A x2\ / a \2 V 2 , 1 , 2 V 2 -V 2 



(A Vl ), = -Av^, ((A, l ) 2 ), = (A, || )^ + -(A, ± ) 2 



ViVj 



(AviAvj)<p - —2 



(A V|| ) 2 - 



(Av ± ) 



i + j (66) 



where v = (v 2 + v 2 + w 2 ) 1 / 2 and 
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Av\\ = 2v 



2„,4 



1 + 



l 2 v 



G 2 {m\ + m2) 2 



-1 




= 2v 





' fa 2 

,^4 



(67) 



Here, again, £ is the impact parameter. 
Using these expressions we find that 



(A(i 2 )>* = ^ 



(Aw,) 2 - 2uAwm 9 9n (Avi) 2 , 9 
(Avi,) 2 - 2vAt> 



i« 2 + 



(A^) 2 
2v 2 



2v 2 

(Vy + vl) 



_ N Awn v 2 + 4v 2 

e- Ae = g- , 

v U 2 rf 



1 H 



(i-Ai) = 



Av|| ti 2 



((A/,) 2 



ft 2 ^ 



(A-|,) 2 ^ | ( A ^) 2 ^ + ^ 



2v„ 



((e-Ae)A^ = — 



2 ^+4j (At; ± ) 2 3^ + 4^ 2 
( At, l|) 13 + o 72 



((i.Ai)A^ = ^ 



21 



(68) 



We integrate these expressions over / in the range from to l max - In doing this we keep only 
the leading terms in the logarithmic factor which appears when we integrate Av\\ and (Afj_) 2 over 
I. In this approximation terms proportional to (Avh) 2 vanish. As a result we find that 



{A(e 2 )) T<u} = C(h, e,i) 2e 2 + hi 2 - (15/4)/j 2 
(e • Ae) T>aJ = -C(h,e,i)e 2 
(A(i 2 )) r , w = C(h,e,i) [e 2 - 2i 2 - (3/4)h 2 

{i-Ai) TtU1 = -C(h,e,i)i 2 , 
(A~h) T ,u, = -C(h,e,i)h, 
((Ah) 2 ) TtUJ = 4C(h,e;i) (e 2 + i 2 -h 2 ), 

((e • Ae)Ah) TiU1 = C{h, e, ~i)h (?,e 2 + 4? - 
<(i- Ai) Ah) TtU , = -C(h,e,i)hi 2 , 



3h 2 ), 



(69) 
(70) 
(71) 

(72) 
(73) 
(74) 

(75) 
(76) 
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where 

C(h, e, i) 



In A 



3/2 



1 + 



In A 



(77) 



i\h\Ve 2 -h 2 e 2 + i 2 -(3/4)h 2 
and Coulomb factor 

A = I ma xv 2 /G(m 1 +m 2 ). (78) 

The upper cutoff l max is determined by the disk dimensions and we set it to be of the order of the 
disk thickness irn (Stewart & Wetherill 1989; SI). Substituting l max = irn and (60) into equation 
(78) we find that 

A = i( Cl e 2 + c 2 i 2 ) > 1, (79) 

where are some constants. By introducing these constants we avoid the need of averaging the 
logarithmic factors over e and i; instead we fix them using numerical data [see Ohtsuki et al. (2002) 
for a similar treatment of Coulomb logarithms] . One can see that all the scattering coefficients that 
we retain in equations (49)-(51) and (52)-(55) are proportional to In A ^ 1. The coefficients of the 
third and higher orders in the Fokker-Planck expansion do not contain this multiplier. Thus our 
neglect of these higher-order coefficients has a relative accuracy of O ((In A) -1 ). 

Using equations (69)-(76) one can easily check that 

<A(e 2 ) + A(i 2 )-^A(P)) T , w = 0, (80) 

<(e • Ae)Ah + (i • Ai)Ah - ^h(Ah) 2 ) T>UJ = 0. (81) 

The first equation is implied by the conservation of Jacobi constant (15). The second equation is 
a statement of {AJAh) T>U) = up to the second order in perturbed quantities [or up to the factors 
~ (In A)" 1 ]. 

Now we perform the final step in our programme — we average the scattering coefficients 
over the Gaussian distribution of relative eccentricity and inclination given by (57). To do this we 
have to perform a series of straightforward but rather lengthy steps which is described in detail 
in Appendix B. Coulomb logarithm In A is always treated as a constant, which is justified by its 
weak dependence on e and i. Of course one can do better than this (see Ida et al. 1993) but for 
our present purposes such an accuracy is sufficient. 

As a result we obtain the following formulae for our scattering coefficients: 

(A(e 2 )) = K{h,d e ,di) (~ttf + ^0 " \ul) , (82) 

(e • Ae) = -K(h,a e ,ai) (V ° + + ^i) - (83) 

(A(i 2 )) = K(h,a e ,^ (\U° - i^ 1 + , (84) 
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(i-Ai) = -K(h,a e ,a i ) (|^o ~ » 



(85) 



{Ah) = - K{kd y^ Ul (86) 
h 

((Ah) 2 ) = K(h,a e ,a i )u£, (87) 



(e 2 Ah) = -hK(h,a e ,ai) (V ° + + ±U^J , 
{i 2 A~h) = -hK(h,a e ,ai) (±U£ - ^uA , 
(e 2 (Ah) 2 ) = ~h 2 K(h,a e ,ai) (u£ + l -U 2 + ^U 2 ^j , 
(i 2 (Ah) 2 ) = h 2 K(h,a e ,a t ) (\u 2 - l -U 2 \ , (91) 



(88) 
(89) 
(90) 



<(e • Ae)Ah) = hK(h, a e , a t ) Q^ 1 - ^uA , (92) 



<(i • Ai)A~h) = -hK(h, a e , at) Qf/o 1 - ±U^J . 
In these expressions we use the following notation: 



(93) 



K{hrWi) = \^e-*l**\ 04) 

a e cr i 



oo 



V} K,a,) = Jdt J^^^ {al+al)t I P [I (<*i " «e) t 




(95) 



(I p is a modified Bessel function of order p) and it is always assumed in equations (82)- (92) that 

Coefficients ((Ae) 2 ) and ((Ai) 2 ) can be trivially computed from (A(e 2 )), (e • Ae), (A(i 2 )), and 
(i • Ai). We also use for the factor inside the Coulomb logarithm the following expression: A = 
Oi(c\d 2 e + c 2 of) > 1. 

We perform some checks and comparisons of our scattering coefficients with the results obtained 
by other authors. First of all, one can check again in a manner analogous to equations (80) and 
(81) that (82)-(93) preserve the Jacobi constant. Second we have checked that our formulae for 
(A(e 2 )) rA ,, (e • Ae) T>w , (A(i 2 )) T ^, (i • Ai) T)W and (Ah) TjUI agree with SI and Ida et al. (2000). Third, 
in the case of a homogeneous planetesimal disk integrals of (A(e 2 )), (A(i 2 )), (e • Ae), and (i • Ai) 
over (3/2)\h\dh should reproduce the averaged viscous stirring and dynamical friction coefficients 
{Pvs}, (Qvs) and {Pdf)i {Qdf) of SI. We numerically checked that this is really the case. 

This completes our calculation of the contribution of Region 1 of e — i space (e > \h\, arbitrary 
i) to the scattering coefficients in the dispersion-dominated regime. 
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5. Discussion and summary. 

We have derived a self-consistent set of equations describing the coupled evolution of the surface 
density and kinematic properties of a planetesimal disk driven by gravitational encounters between 
planetesimals. The assumption of a dispersion-dominated velocity regime is used throughout the 
calculation, which is reasonable for planetesimal disks in their late evolutionary stages. Thus this 
paper serves as a logical extension of Paper I which was devoted to the study of the shear-dominated 
velocity regime. 

The evolution equations (49) and (52) are of advection-diffusion type; the coefficients entering 
them are nonlocal which is in contrast with previous results (Paper I; Ohtsuki Sz Tanaka 2002). 
This is a natural outcome of the scattering in the dispersion-dominated regime since planetesimals 
can explore large portions of the disk (compared to their Hill radii) in the course of their epicyclic 
motion. 

We have also derived analytical expressions for the scattering coefficients entering the different 
terms of evolution equations (50), (51), (53)-(55). The analytical treatment was enabled by the 
use of the two-body scattering approximation which becomes valid in the dispersion-dominated 
regime. Our expressions are accurate up to fractional errors ~ (InA) -1 , where In A 3> 1 is a 
Coulomb logarithm. Following the methods developed in SI and Ohtsuki et al. (2002) for the 
case of homogeneous planetesimal disks it might be possible to improve our calculations by taking 
subdominant higher order terms into account (they contribute typically at the level ~ 10% — 20% 
and are neglected in our present consideration). 

Using this system of evolution equations supplemented with the information about the be- 
havior of the scattering coefficients one can self-consistently study the evolution of inhomogeneous 
planetesimal disks. Arbitrary distribution of masses of interacting planetesimals is allowed for but 
the evolution of mass spectrum is not considered in the present study. Thus our equations describe 
the disk evolution on timescales short compared with the timescale of the mass spectrum evolution. 
This should be a good assumption for studying effects on the disk caused by the gravity of massive 
protoplanetary embryos (because they can induce changes of the disk properties on rather short 
timescales). This deficiency can also be easily remedied in the future when the need to study very 
long-term disk evolution arises. Our approach can naturally incorporate physical mechanisms other 
than just gravitational scattering, for example gas drag or migration [see Tanaka & Ida (1999)] and 
we are going to study their effects in the future. 

In the following paper (Paper II) we describe the embryo-planetesimal scattering and derive 
equations governing this process in various velocity regimes. This would allow us to provide a 
complete description of the disk evolution caused by both the presence of isolated massive bodies 
and the continuous sea of planetesimals. 

I am indebted to my advisor, Scott Tremaine, for his patient guidance and valuable advice 
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which I received during the work on this paper. Financial support provided by the Charlotte 
Elizabeth Procter Fellowship and NASA grant NAG5-10456 is gratefully acknowledged. 
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A. Calculation of the instantaneous surface density. 

Suppose that we know the behavior of the surface density of planetesimal guiding centers N(h) 
in the whole disk as well as the random velocity distribution function ip(e, i, h) (we will initially 
carry out calculations for the case of general velocity distribution function). We want to compute 
the instantaneous surface density N mst (xo) at some point xq, where xo is scaled by the reference 
radius ao and is thus dimensionless. 

To perform the desired conversion let us denote the number of planetesimals per unit dv x dv y dv z dxodz 
as g tnst (v x ,v y , v z , xo, z) and the same number per unit dedidhdrdu as g(e, i, h, r, u). The Jacobian 
of the transformation between these sets of coordinates [see equations (2) and (13)] is 



J = ^tfei 
2 

(z and velocities v x ,v y ,v z are assumed to be dimensional). Clearly, 

2g(e,i,h,T,u) 



N inst (x ) = J g inst (v x ,Vy,v z ,x ,z)dv x dv y dv z dz = J 



aiQ 3 ei 



dv x dv y dv z dz. 



(Al) 



(A2) 



To compute N mst one needs the following expressions obtained using (2) and (13): 



„ vl + Avl .„ v 2 + tt 2 z 2 , 2v v 
e 2 = x — J , i 2 = -S— — — , h = x + y 



(n ao ) 2 ' 



(ftao) 2 ' 



(A3) 



Now we specify a distribution function (17) and find that 



N mst (x ) 



2vr 2 a^ 3 



n dv x dv v dv z dz 
N(h) Z" 2 " 
a 2 e{h)o-t{h) 



cxp 



1 (vl + Av 2 v 2 + n 2 z 2 



n 2 a 2 \ 2a 2 



2af 



(A4) 



with h given by (A3). We can easily perform the integration over dv x dv z dz. Switching also from 
v y to h we finally obtain that 



oo 

V2TT J (Teih) 



(x - h) 



21 



2ol 



(A5) 



Note that N 



inst j 



xo) 



N(x ) when a e 0. Also, if N(h) = const, then N inst (x ) = N 



according to formula (A5). 



B. Averaging over e, i distribution. 



We illustrate the procedure of averaging the scattering coefficients over the distribution func- 
tion (57), using the coefficient (e • Ae) as an example. Calculation of other scattering coefficients 
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is analogous. Using equations (70) and (94) we write: 

8 In A 7 7e 3 dediexp [-e 2 /(2a 2 ) - i 2 /{2d 2 )] 



(e-Ae) = -- 



3ir\h\a 2 a 2 J J JpZtf \~ e 2 + -2 _ (3/4)^2 



3/2 ' 



(Bl) 



The lower limit of the integration over e is set to h because only test bodies with eccentricities 
bigger than this value can experience close encounters with the reference particle. Introducing new 
variables x and y by 



h 2 ^ ' /l 2 



we can rewrite the double integral in (Bl) as 



00 00 







1 + (l/4)x 2 

(l +x 2 +y 2)3/2' 



-a|a; 2 — a?i/ 2 



dxdy, 



where 



2 _1^ 2_i^ 



Making another change of variables x = y^cos^, y = y/rsm.0 we obtain that 



00 



(l +r )3/2 



exp 



tt/2 



r r cos 26* 
1 + - H I exp 



r cos 20 



(a? " 



d(9. 



(B2) 



(B3) 



(B4) 



(B5) 



Using the identity (Gradshteyn & Ryzhik 1980) 

J e zcost cos nt dt = irl n (z), (B6) 


where I n is a modified Bessel function, and introducing definitions (95) and (96) we finally obtain 

(B7) 



\h\e-*/&i) ( £7° + Itf + \U\ ) . 



71", r, _l2 //0 -2 
2 



Using this result and (Bl) we finally arrive at equation (83). 



